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ABSTRACT 

Presented in this paper is a derivation of a 2D catalytic reaction-based model to solve combinatorial 
optimization problems (COPs). The simulated catalytic reactions, a computational metaphor, occurs 
in an artificial chemical reactor that finds near-optimal solutions to COPs. The artificial environment 
is governed by catalytic reactions that can alter the structure of artificial molecular elements. Altering 
the molecular structure means finding new solutions to the COP. The molecular mass of the elements 
was considered as a measure of goodness of fit of the solutions. Several data structures and matrices 
were used to record the directions and locations of the molecules. These provided the model the 
2D topology. The Traveling Salesperson Problem (TSP) was used as a working example. The 
performance of the model in finding a solution for the TSP was compared to the performance of 
a topology-less model. Experimental results show that the 2D model performs better than the 
topology-less one. 
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1. Introduction 

Solutions to combinatorial optimization problems (COPs) have practical real-world importance 
because most real-world problems are combinatorial in nature. Most COPs have been shown to 
be AAP-complete. Exact algorithms have been proposed to these problems but prove inefficient 
for large problem instances [14]. Craph-based algorithms such as branch and bound [24], as 
well as distributed multi-agent based heuristics such as genetic algorithms [22], memetic algo¬ 
rithms [20, 10, 11], tabu search [25], simulated annealing [18], simulated jumping [2], neural 
networks [19], and swarm intelligence [12, 13, 8] have been used to find time-restrained optimal 
and near optimal solutions for these problems. 

In recent years, the chemical systems of living organisms have been shown to possess inherent 
computational properties [15, 1, 3]. This discovery provided researchers the chemical metaphor 
as a paradigm for computation [6, 9, 4, 16, 5, 7]. Under this computational framework, molecules 
are considered as solutions, while interactions among molecules represent computational proce¬ 
dures. 

Presented in this paper are the mapping of permutations to molecules, and the derivation 
of two stochastic functions that model catalytic reactions. The two functions simulate a unary 
catalytic reaction and a binary catalytic reaction. These reactions create new molecules where 
the average goodness of fit of the product is better than the average goodness of fit of the 
reactants. A molecule encodes a permutation such that the unary catalytic reaction reorders 
the element in the permutation while the binary catalytic reaction creates new permutations. 
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2. Model Development 

This section discusses the development of the 2D artificial catalytic reactor (2DACR). The 
2DACR models the dynamics of artificial molecules in a 2D environment. The environment 
is driven by several rules of interactions to produce a set of optimal or near-optimal solutions 
to a COP. The 2DACR is defined by a triple 2DACR(M, R, A), where M is a set of artificial 
molecules, R is a set of reaction rules describing the interaction among molecules, and A is 
an algorithm driving the reactor. In this paper, the molecules in M are permutations while 
the rules in R are reordering algorithms that create new molecules. The algorithm A describes 
how the rules are applied to a vessel of artificial molecules simulating a well-stirred, 2D reactor. 
Further, the 2DACR is partitioned by A into different levels of reaction activities. The level of 
reaction activity is a function of molecular mass. 

2.1. Mapping 2DACR to COP 

In a particular COP, if n is the length of a permutation, then all n-sized permutations n„ 
are molecules. The inherent value of a permutation vr G n„ is the mass of the molecule and 
is computed by the objective function of the COP. For example, if the COP is a traveling 
salesperson problem (TSP) with a set of cities V, a set of paths R,and a cost matrix G, then 
all permutations of the |P| cities are molecules. Each n-sized molecule is an encoding of a 
Hamiltonian tour such that the ordering of the cities Uj G P represents a molecular structure. 
Each city Vi is a distinct atom in the molecule. The cost of traversing a Hamiltonian tour 
is the mass of the molecule and is a measure of goodness of fit of the Hamiltonian tour. If the 
objective function is a minimization of the tour cost, then light molecules encode the desirable 
solutions. In this example, is defined in Eq. (2.1) where gij £ G is the cost measure associated 
with path {vi,Vj) G E: 

n—1 

fv — 9n,l + 'y ( 2 - 1 ) 

i=l 

2.2. Artificial Catalytic Reactions 

If two molecules mi and m 2 collide, they react following a catalytic reaction of the form 
mi -|- m2 -|- C —mi -|- m2 -|- m3 -|- m 4 , 

where m 3 and m 4 are product molecules and C is a catalyst. The reaction is a mathematical 
function 

Ri : M X M —M X M X M X M, 

where M = {mi Vi = 1, 2,... , |n„|}, m* is a molecule and is a linear data structure representing 
a permutation, and |n„| is the cardinality of the solution space of the COP (i.e., \M\ = |n„|). Ri 
performs a reordering of solutions and is dependent on an n-long vector S with binary elements. 
The elements s £ S are computed following this algorithm; 

1. Let random(n) be a function that returns a random integer between 1 and n. 

2. Let index(mi, m 2 j) be a function that returns the position of atom m 2 j' in molecule mi. 

3. Let oji be a set of distinct atoms in molecule mi. 

4. Let mij be the jth atom in molecule mi. 

5. Let the integers i and j be indeces 

6 . Initialize the vector S = 0. 
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7. Initialize ui = {}. 

8 . Set i = 0. 

9. Let j = random(n) be the point of collision between mi and m 2 . 

10. Do the following until i = n or m 2 j G coi- 

(a) Increment i by one. 

(b) Set oji = uji 

(c) Set Sj = 1. 

(d) If m 2 j ^ oJi, then j = index(mi, m 2 j). 

Once S is obtained, the molecules m 3 and can be computed using the following algorithm: 

1. For each i = 1,..., n do: 

(a) If Sj = 1, then set ma^j = mi^j and myi = m 2 ,i; 

(b) Else set set m-^^i = m 2 ,i and = mi^*. 

If a molecule ms hits the bottom or walls of the reactor, a catalytic reaction of the form 
ms + C —mQ 

happens. The reaction is a mathematical function 
i?2 : M —> M 

described by the following algorithm: 

1. Let j = random(n) be the point of collision between ms and the reactor wall or bottom. 

2. Let i = random(2) — 1. 

3. If i = 0, then do the following: 

(a) If j = n, then me,! = 

(b) Else m 6 ,j+i = 

4. If i = 1, then do the following: 

(a) If j = I, then m 6 ,n = 

(b) Else mQj-i = ms,j. 

2.3. Two-Dimensional Reactor 

The reactor algorithm A operates on a matrix T of molecules and catalysts, where T = = 

V = 0 Vi = 1,..., / A j = I,..., J}, / > 1, J > 1 and [I x J) << \M\. mu is the /cth 
molecule while the catalyst is when tij = 0. Here, the matrix element tij serves a placeholder 
for the molecule such that the expression tij = should be understood as “m^ is at Uj.'" 

2 . 4 . Molecule Direction 

Associated with each molecule m^ is a direction of the molecule. The range of values for 
is dependent on the stencil used by A. A stencil is a set of directions from an element in a grid 
environment. In a 2D environment, the possible stencils are 5-point stencil and 9-point stencil 
(Figure I). There are five possible directions in a 5-point stencil: no movement (0), due North 
(N), due East (E), due West (W), and due South (S). In a 9-point stencil, the nine possible 
directions are 0, N, NE, E, SE, S, SW, W, and NW. An arbitrary integer may be assigned to 
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Figure 1: Two of the possible direction stencils for a 2D environment. 


each of these directions reserving the value zero for the no movement direction. Elements at the 
borders and at the corners also follow any of these stencils. 

A random direction is assigned to each molecule during the start of the simulation. After a 
collision of type Ri and following a 9-point stencil, the directions of the products are determined 
as follows; 

1. If the mass of a molecule is greater than the average mass of the products, then the direction 
assigned is random(NW, N, NE). 

2. If the mass of a molecule is lesser than the average mass of the products, then the direction 
assigned is random(SW, S, SE). 

3. If the mass of a molecule is the same^ as the average mass of the products, then the 
direction assigned is random(W, 0, E). 

After a collision of type R 2 and following a 9-point stencil, the direction of the product is 
determined as follows: 

1. If the mass of the product is greater than the mass of the reactant, 
assigned is random(NW,N,NE). 

2. If the mass of the product is greater than the mass of the reactant, 
assigned is random(SW, S, SE). 

3. If the mass of the product is the same as the mass of the reactant, 
assigned is random(W, 0, E). 

2.5. Collision Matrix 

The algorithm A also maintains a collision matrix C whose elements are defined by Eq. (2.2). 
The matrix C has the same dimension as T and records which molecules will collide at which 
element in T. Two molecules mk and mi will collide at tij if ctj = {mj,mk). The collision 
will obey the reaction rule defined by i?i. A molecule will collide with the border at tij if 
Cij = mk- In this case, the reaction rule R^ will be applied. If Cjj = 0, then no collision will 
occur at Uj. Each element Cij G C is updated by tracing the movement of a molecule via its 
direction. Usually, C is updated via any of the two popular computational order; row-major 
ordering and column-major ordering. However, these two popular ordering techniques have 
inherent biases. The row-major ordering scheme is biased towards molecules at lower values of 
i (i.e., at the upper portion of T) while the column-major ordering is biased towards those at 
lower values of j (i.e., at the left portion of T). To remove the biases brought about by these 
two ordering schemes, ordering schemes based on space-filling curves are recommended. In this 


then the direction 

then the direction 

then the direction 


^In practical application, this may mean as statistically the same. 
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Figure 2: Order of computation for a 4 x 4 T following the Morton ordering scheme. 


paper, the Morton ordering scheme is used as a working example and is discussed in the next 
section. 


{ {mk,mi) if rrik A m/ collide at Uj 

mk if rrik collides with the wall at tij (2-2) 

0 if no collision occurs at tij 

2.6. Ordering of Computation 

The Morton ordering (Figure 2) is an ordering scheme based on Peano-Hilbert plane-filling 
curves. A plane-filling curve is a curve drawn on a plane and fills it. A plane-filling curve is 
efficiently generated using a recursive algorithm that divides the plane at each recursion level. 
The principle followed at each recursion is the self-similarity principle such that the structure 
of the curve at the higher level is the same as the structure of the curve at the lower level. 
Traversing a curve means to enumerate the points along the curve. When the curve fills the 
plane, the traversal of the curve also means the traversal of the plane at the order defined by 
the curve. The order of traversal is defined by a universal turtle traversal algorithm [17], which 
traverses a self-similar space-filling curve based on a movement specification table. 

2. 7. Evolution of the Reactor 

The evolution in A is realized by applying the following algorithm: 

1. Initialize T with I x J molecules selected randomly from M. 

2. Compute the molecular mass of each mk G T and initialize each with a direction. 

3. Compute for C using the Morton ordering scheme. 

4. Apply the reaction rule Ri for any mi and m 2 . 

5. Apply the reaction rule R 2 for heavy molecules that collide with the reactor walls and 
bottom. 

6. Decay the heavier molecules by removing them out of T and replacing them with randomly 
selected molecules from M. 

7. Repeat steps 2 to 5 until T is saturated with lighter molecules. 
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One iteration of A constitutes one epoch in the artificial reactor. The sampling procedure 
gives molecules with low molecular mass a higher probability to react or collide with other 
molecules. This mimics the level of excitation energy the molecule needs to overcome for it to 
react with another molecule. This means that the lighter the molecule, the higher the chance 
that it will collide with other molecules. Step 6 of algorithm A requires a metric for measuring 
saturation of molecules. The reactor is considered saturated when T has no more 0 element 
(i.e., the catalyst is already exhausted). 


3. Experimental Results 

A 2DACR was run to solve an instance of a symmetric TSP. The 2DACR used the 5-point 
stencil and the Morton ordering scheme. To assess the performance of the 2DACR, a topology¬ 
less artificial chemical reactor (ODACR) recently employed by other researchers [23, 21] was 
also run to solve the same TSP instance. A single-processor Pentium IV machine with 1.2GHz 
bus speed running under a multiprogramming operating system was used to run the 2DACR 
and the ODACR simulations. The simulations were repeated 10 times while the best minimum 
Hamiltonian tour length for each run were recorded. The values recorded were averaged and 
the standard deviation computed. Shorter tour lengths imply better tour costs and are much 
desirable. To remove the initialization bias, both 2DACR and ODACR used the same initial set 
of 500 molecules, with the exception that the molecules in 2DACR have their respective initial 
locations and directions while those of the ODACR have none. The 2DACR utilized a 30 x 30 
matrix T with 400 elements in T acting as the catalyst C. Table 1 compares the average tour 
lengths found by 2DACR and ODACR on five sets of random instances of symmetric 50-city 
TSPs. The table shows the average value of 10 runs for both 2DACR and ODACR with their 
respective standard deviations. Based on the result presented, it can be seen that 2DACR’s 
performance is better than the performance of the ODACR employed by other researchers. 


Table 1: Comparison of average Hamiltonian tour length found by 2DACR and ODACR on five 
sets of random instances of symmetric 50-city TSPs. The values are averaged over 10 runs. The 
values in parenthesis are the standard deviations. 


Problem 

ODACR (std. dev. 

) 2DACR 

1 

5.89 (0.40) 

5.87 (0.45) 

2 

6.17 (0.08) 

6.15 (0.11) 

3 

5.65 (0.21) 

5.59 (0.18) 

4 

5.70 (0.98) 

5.67 (0.77) 

5 

6.15 (0.54) 

6.14 (0.55) 


4. Concluding Remarks 

An algorithm that models catalytic reactions in 2D was designed to solve COPs using the TSP 
as a working example. Solutions to an instance of TSP via 2DACR were found better on the 
average than those found by ODACR. Molecular directions and locations were incorporated that 
provide the model with a 2D topology. The order of computation via Morton ordering might have 
removed the bias inherent in row-major and column-major ordering schemes. However, more 
work is needed that will compare the performance of 2DACR using these ordering schemes. 
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